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ABSTRACT 

The  application  of  Mitrovic's  Method  enables  the  control  systems 
engineer  to  analyze  and  synthesize  control  systems  as  a  function  of  two 
parameters.   The  feedback  compensation  of  sampled-data  control  systems 
is  thoroughly  investigated  as  a  multi-parameter  problem,  with  emphasis 
placed  on  the  design  of  the  compensation  in  the  Mitrovic  parameter  plane. 
The  minimization  of  the  control  system  settling  time  is  shown  to  be  parti- 
cularly adapted  to  the  parameter  plane.   Several  examples  are  included 
which  demonstrates  this  usefulness  of  the  parameter  plane  in  solving  the 
minimum  settling  time  problem.  A  stability  relationship  is  devised  for 
digital  computer  determination  of  dynamic  system  stability  which  permits 
the  formulation  of  an  n-dimensioned  parameter  space.   Applications  of  this 
digital  creation  of  an  n-dimensioned  space  are  considered.   This  ability 
to  design  feedback  compensation  for  the  sampled-data  system  as  a  function 
of  more  than  one  parameter  is  shown  to  lead  to  greater  accuracy  in  closed 
loop  root  positioning;  and  in  certain  cases,  to  permit  the  design  con- 
sideration of  several  constraints  simultaneously. 
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CHAPTER  T 
GENERAL  CONCEPTS  AND  THEORY 
1.1   Introduction, 

The  design  of  multivariable  feedback  compensation  for  a  feedback 
control  system,  linear  or  sampled-data,  has  heretofore  been  accomplished 
largely  by  one  of  two  methods.   Either  repeated  calculations  using  the 
classical  one  parameter  design  techniques,  such  as  root  locus,  has  been 
used;  or,  if  a  more  accurate  study  is  required,  the  designer  could  resort 
to  analog  simulation  and  a  "trial  and  error"  approach  to  the  design  pro- 
blem.  Dusan  Mitrovic,  [3],  in  his  early  work  with  the  algebraic  methods 
for  locating  the  roots  of  a  polynomial  introduced  a  new  analysis  technique 
that  permitted  analysis  of  the  root  location  in  terms  of  two  variable  co- 
efficients.  Later,  Siljak,  [4]  and  [5],  devised  additional  algebraic 
methods  whereby  two  parameters  could  be  selected  from  a  portion  of  any, 
cr  all,  of  the  polynomial  coefficients,  and  the  root  analysis  performed 
in  terms  of  these  rarameters.   This  extension  has  proven  exceptionally 
useful  in  two  parameter  analysis  of  control  systems  whose  characteristic 
equations  have  coefficients  that  are  each  functions  of  the  same  two  vari- 
ables or  system  parameters. 

The  object  of  this  thesis,  then  is  to  apply  the  results  of  Mitrovic 
and  Siljak  to  the  design  of  feedback  compensation  for  the  sampled-data 
control  system.   A  secondary  objective  is  to  show  the  superiority  of  the 
two  parameter  method  over  the  more  classical  techniques,  and  to  demon- 
strate the  adaptability  of  the  entire  design  problem  to  solution  by  digital 
computer. 


i  .  2   8  -!  s  i  c  Equ  a  t  j.  on  s . 

Let  us  consider  the  characteristic  equation  of  the  sampled-data  feed- 
back control  system  as, 

F(z)  -  1  +  GH(z)  =  0  (la) 

or,  F(z)  =  2_  akz  =  0         ,  (lb) 

KrO 

where  the  a  '  s  are  constant  or  parameter  varying  coefficients,  and  z  is 

one  of  the  complex  z-plane  locations  of  the  closed  loop  roots  of  the  system. 

Now,  consider  the  generalized  point  in  the  z-plane, 

whe  re  |  z  j  =    W£ 

and  arg   z  =  e    r 

■  cos  frf=j   sin  J# 

-  £  +  «•  -  r;  >* 

If  this  general  coordinate  is  substituted  into  equation  (lb),  the  character- 
istic equation  may  then  be  written  as: 


K=0  *    J  . 


(3) 


Equation  (3)  is,  in  general,  a  non-zero  quantity  and  assumes  a  zero  value 

only  when  the  value  of  6u?  and  C  assumed  corresponds  to  that  of  one  of 

the  closed  loop  roots  of  the  system.   Thus,  if  we  require  that  F(z)  =  0, 
equation  (3)  describes  the  closed  loop  pole  locations  if  the  a,  's  are 

constants,  and  every  point  on  the  z-plane  root  locus  if  one  of  the  a,  's 
is  parameter  varying.   Therefore, 

.2  vfc  v  k   ~ 


£\<<o*£    +  J4yi-Sfr> 


)  (4) 

K--o 
mans  the  characteristic  equation  from  the  z-plane  to  the  Mitrovic  plane  if 

the  a  's  are  constants,  and  the  parameter  plane  if  the  a,  's  are  variable. 
k  K 

2 


■ 
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Only  the  parameter  plane  will  be  discussed  hereafter. 

i  .  3  Con  form  a  I  '  ap  -in^  Techniques . 

The  requirement  that  equation  (4)  be  identically  equal  to  zero 
implies  that  both  the  real  and  imaginary  parts  are  zero.   Siljak,  [5], 
demonstrated  that  this  magnitude  requirement  on  equation  (4)  leads  to 
the  two  simultaneous  equations, 

tl  ak  ^  V  £.  >  -  °  <5a> 

and  , 

£  ax  u)"  v  Si >  *  °     *  (5b) 

where  T,  (   C  )  nnd  U  (  ^  )  are  the  Chebyshev  functions  of  the  first  and 
second  kind  respectively.   Furthermore,  since 


T 
k 


:<<^>   =      ^Uk<>,>    "  Vl(£>  '  (6) 

equations    (5a)    and    (5b)    can  be   rewritten   as 


2_-  s4  Vi<^>  -°  <7a> 

and 


respectively,  with  the  U,  (  v^  )  s  generated  by  the  recursive  Chebyshev 
formula 


(7b) 


Uk+1 


(SL>-2£  v^  -uk-i(Q>     <8> 


with  UQ(  ^  )  ^  0 

ut(  Cj        =  l 
*  -l 


u-i<^> 

Siljak,  [4]  and  [5],  considered  the  case  in  which  the  a,  '  s  are  linear 
functions  of  two  scalar  parameters, 

3 


ak=o<bk+^ck  +  dk  .  (9) 

The   substitution  of  equation   (9)    into  equations    (7a)    and   (7b)    forms   the 
set  of   simultaneous   equations, 


and 


where 


^V^S^/^V^^V^S^0  • 


(10b) 


<W\4>v 


Bi=   2_  -bk4  Dk-i<£>    •    B2-   z_   bk  4  \<S,' 

s-  fT -ck  4  °k.i<^>   •  c2=  2T  ck^Kvse) 

Equations  (10a)  and  (10b)  can  then  be  solved  simultaneously  for  the  two 
parameters,  alpha  and  beta  as; 

C1D2  "  C2D1               •  Q             B2D1  "  B1D2      .        (11) 
^s  ft-   

B1C2  '  B2C1  B1C2  "  B2C1 

Equations  (11)  are  the  fundamental  results  of  Siljak,  [4]  and  [5], 

and  provide  the  basis  for  the  work  which  follows.   Equation  (11)  can  be 
interpreted  as  the  calculated  values  of  the  two  parameters  alpha  and  beta 
required  to  put  a  closed  loop  root  at  a  particularly  chosen  location  on 
the  z-plane.   Thus,  by  fixing  6JL  and  varying  S~   from  -1  to  +1,  upper 
semi-circle  contours  of  constant  6j^  can  be  mapped  on  the  parameter  plane. 
Thus,  the  argument  of  each  upper  semi-circle  complex  root  conformally  maps 
onto  a  single  point  on  the  parameter  plane,  the  coordinates  of  which  are 
the  values  of  the  two  parameters  producing  that  particular  root  configura- 
tion.  Similiarly,  in  equation  (11),  ^    can  be  fixed  while  C^^is   varied 
from  zero  to  infinity  producing  loci  of  constant  O  ,  or  angle,  on  the 
parameter  plane. 


.     ■•• 


• 


A  third  type  of  conformal  mapping  is  yet  required  to  completely  specify 
all  the  roots  of  the  characteristic  equation  on  the  parameter  plane.   The 
loci  of  constant  Cj     and  constant  ^  determine  all  the  complex  roots. 
The  loci  of  points  on  the  parameter  plane  corresponding  to  constant  real 
roots  on  the  z-plane  can  be  obtained  directly  from  the  characteristic  equa- 
tion, equation  (lb).   The  substitution  of  equation  (9)  into  equation  (lb), 

and  letting  z  =  C^  ,  yields 

P  Q 


o<  ZbkCrK  +  /3j\o7K  +  5Ldk<r;K    -o 


(12) 

**°  'k^j  K--0 


which  for  a  given  OZ  is  a  straight  line  on  the  parameter  plane,  regard 
less  of  the  order  of  the  characteristic  equation. 
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CHAPTER  II 
SECOND  ORDER  SYSTEMS 
The  second  order  sampled-data  control  system  provides  an  excellent 
model  for  developing  the  feedback  compensation  techniques  from  applica- 
tions of  the  equations  of  Chapter  I.   For  the  initial  design  application, 
consider  the  pure  second  order  system  shown  in  figure  (la).   The  selected 
design  goal  is  to  minimize  the  settling  time  by  suitable  choice  of  the  feed- 
back parameters,  a.,  and  a».   If  the  states  of  the  system  of  figure  (la^  are 


defined  as: 


X^t)  =  c(t) 


and  X2(t)  5  c(t)        ,  (13) 

then  the  system  can  be  represented  by  the  signal  flow-graph  of  figure  (lb), 
and  the  vector-matrix  equation 

X[(k  +  1)T]  =  ^   X(kT)  +  A  R(kT), 


(14) 


where 


I    = 


(a  KTT)/2 


(a^T) 


T 

1 


(a2Rf)/2 


(a2KT) 


and 


A  = 


KT2/2 


KT 


(15) 


(16) 


The  characteristic  equation  of  the  system  is; 


2z2  +  z(2a  KT  +  a  KT2  -  4)  +  (2  -  2a  KT  +  a^T2)  =  0 


2.1  Constant  Magnitude  Loci. 

The  duration  of  the  transient  response,  or  settling  time  for  a 
linear  system  of  the  second  order  is,  [6], 


• 


<&§                 V' 

Figure  1                    ' 
Pure   Second  Order  System 
Diagrams 
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(b)    Signal  Plow-Graph 
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For  a  stable  system  then,  the  contours  of  lines  of  constant  settling  time 
of  a  second  order  system  in  the  s-plane  are 

or  s  =  -  4/t  .  (18) 

s 

Accordingly,  from  one  definition  of  the  z-transformation,  [2], 

z  =  eTs  ,  (19) 

we  can  determine  corresponding  contours  of  constant  settling  time  in  the 

z-plane  as, 

z   =  e"4T/ts  ,  (20) 

which,  since  T  and  t  are  constants  for  any  particular  design,  are  circles 

s 

-4T/t 
centered  at  the  origin  of  the  z-plane  of  radius  e     s. 

Thus,  the  loci  of  constant  60^  contours  on  the  parameter  plane  can 

be  interpreted  as  loci  of  constant  settling  times;  and,  when  C^z   and  T 

are  specified,  the  appropriate  t  can  be  uniquely  determined.   Solving 

5 

equation  (20)  for  t  , 

s 

t  =  -  4T/ln  6JZ         .  (21) 

Equation  (21),  however,  must  be  interpreted  in  terms  of  the  z-plane,  dis- 
crete, characteristics.   For  example,  for  a  step  input,  the  minimum  set- 
tling time  obtainable  is  nT  seconds,  for  an  n-th  order  system,  when  all 
the  closed  loop  roots  are  at  the  origin  of  the  z-plane.   Accordingly, 
when  6J2  approaches  zero  in  equation  (21),  t  should  approach  nT.   Thus, 
for  proper  interpretation,  equation  (21)  must  be  written  as 

t  =  nT  -  4T/ln  oX,      .  (22) 

s 

The  particular  conclusion  of  interest  here  is,  in  terms  of  the  parameter 
plane,  to  minimize  t  one  has  only  to  minimize  CJZ  within  the  stable 


' 


region  of  the  plane. 

To  demonstrate  the  effectiveness  of  the  constant  C^-    loci  in  aid- 
ing the  designer  meet  settling  time  criteria,  consider  again  the  pure 
second  order  system  of  figure  (1),  whose  characteristic  equation  is  given 
on  page  6.   It  is  first  necessary  to  define  the  system  parameters,  so  let 

(X   -  axKT2  (23) 

and  "  fi     &   a2KT 

By  reference  to  equation  (9)  and  the  characteristic  equation,  then 

"o  -  2  bo  =  >  co  •  "2 

d  •  -4  b  -  1  c  -  2 

d2  -  2  b2  -  0  c2  =  0 

which  when  substituted  into  equation  (10)  yields 

Bl  =  "U-l<  5.)  B2=  ^Ul<  S  > 

Cl   =   2U-l(S7.  >  C2  =  2^U1<SJ 

Dx   =   -2U_T(  £j    "2   GJjo^  S?) 

D2  -  -4u\uT(  C2)  +  2  tf2  U2(  £) 
These  factors  can  now  be  substituted  into  equation  (11)  and  the  constant 
63  2.     l°ci  plotted  with  the  aid  of  a  desk  calculator  or  a  digital  computer. 
Note  that  if  a  digital  computer  is  used,  as  is  most  likely,  a  limiting 
process  is  required  to  remove  the  common  factor  of  Cu^     from  the  numera- 
tor and  denominator  of  the  parametric  equations  when  attempting  to  locate 
the  O^.  =  0  point  on  the  parameter  plane.   Figure  (2)  shows  the  resulting 
lines  of  constant  C3 g.  °n  the  parameter  plane  for  this  system. 

From  figure  (2),  minimum  settling  time  is  predicted  when  <iV  -  0, 
or  when 

<*-  =  1.0  , 

and  /&      -  1.5 

9 


•         : 
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In  order  to  determine  the  system  response  as  predicted  by  equation  (14), 
Let 

T  =  K  =  1.0  , 

then  a.    =  1.0 

and  a«  =  1.5  is  required. 

Tnble  1  lists  the  calculated  magnitude  values  of  the  states  of  the  system 
at  consecutive  sampling  intervals  using  these  values  of  the  two  parameters, 
The  time  response  of  the  states  is  plotted  on  figure  (3) .   A  unit  step  in- 
put is  assumed. 

TABLE  1 
nT  X^t)  X2(t) 

0  0  0 

1  0.5  1.0 

2  1.0  0 

3  1.0  0 

4  1.0  0 

Since  in  this  example,  n  *  2,  these  results  are  in  agreement  with  that 

predicted  by  equation  (22),  and  the  response  is  the  minimum  prototype, 

ripple  free  response. 

In  order  to  demonstrate  the  most  important  design  advantages  in  using 

the  parameter  plane  for  feedback  compensation;  assume  that,  for  some  other 

reason,  one  of  the  feedback  coefficients  is  constant,  say  a.,  is  fixed  at 

2.25.   Let  K  =  T  =  1.0,  also  be  constant.  With  these  constraints  on  the 

design  problem,  the  parameter  plane  affords  an  excellent  potential  for 

setting  the  remaining  degree  of  freedom  to  minimize  the  settling  time, 

in  so  far  as  possible  with  feedback  compensation  alone.   Referring  to 

figure  (2)  for  <=<  ■  2.25,   &)?  is  minimum  at   /.)  =0.5  when  /&    -   1.875. 

From  equation  (22)  , 

10 
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Figure  3 
Response  of  Pure  Second  Order  System 
with  Minimum  Settling  Time  Compensation 
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t  =  2T  -  4T/ln  0.5 


or 


t   -  7.8T 
s 


is  the  nredicted  settling  time.   Table  2  lists  the  state  variable  magni- 
tudes at  the  indicated  sampling  interval  in  response  to  a  unit  step  input 
utilizing  the  above  parameter  values  in  equation  (14).   Figure  (4)  shows 
the  time  response  of  the  state  variables  of  the  system  as  a  function  of 
the  sampling  interval,  or  time. 
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The  results  indicated  in  Table  2  appear  to  discredit  the  validity  of 
equation  (22)  in  predicting  settling  time.   Actually,  equation  (22)  simply 
nrovides  an  engineering  approximation  to  the  settling  time  problem.   Analy- 
sis of  the  results  of  Table  2  indicates  the  output  position  is  within  1.9570 
of  the  final  value  after  8T  seconds,  although  an  actual  zero  velocity  and 
final  position  of  0.4444  was  not  achieved  until  the  41st.  sampling  interval. 
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TABLE  2 

Xj(t) 

x2(t) 

0 

0 

0.5000 

1 . 0000 

0.5000 

-1.0000 

0.3750 

0.7500 

0.5000 

-0.5000 

0.4063 

0.3125 

0.4688 

-0.1875 

0.4297 

0.1094 

0.4531 

-0.0625 

0.4395 

0.0352 

0.4473 

-0.0195 

0.4429 

0.0107 

0.4453 

-0.0059 

0 . 4440 

0.0032 

0.4447 

-0.0017 

0.4443 

0.0009 

0.4445 

-0.0005 

0.4444 

0.0003 

0.4445 

-0.0001 

0 . 4444 

0.0001 

0.4444 

-0.0000 
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Figure  Ij. 
Response  of  Pure  Second  Order  System 
with  Constraints 
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A  similar  test  was  performed  with  a  =  2.25  and  a  =  1.50  for  6J,  =  0.8. 
Equation  (22)  predicts  the  settling  time  as  t  =  20T,  at  which  time  the 
output  had  settled  to  within  3.3%  of  the  final  value.   In  that  example, 
zero  actual  velocity  did  not  occur  until  t  =  99T.   However,  the  settling 
time  prediction  of  equation  (22)  remains  within  generally  accepted  good 
engineering  design  proximity  of  the  desired  value.   It  must  be  borne  in 
mind  however,  that  for  s     different  from  1.0,  the  conventional  definition 
of  system  error  cannot  be  applied.   Fortunately,  in  many  cases,  input/out- 
put scaling  can  be  applied  which  permits  the  position  voltage  feedback  to 
be  set  at  other  than  unity. 

2.2  Constant  Angle  Loci. 

Equation  (11),  the  parametric  equations  describing 

,     C1D2  "  C2D1  >  A  B2D1  "  B1D2  (U> 

B1C2  "  B2C1  B1C2  '  B2C1 

the  loci  of  the  complex  roots  on  the  parameter  plane  are  functions  of  both 

CJz      and   C^  .   The  constant  magnitude  curves  were  obtained  by  fixing 
C02   and  allowing  JT  to  vary  as  the  mapping  parameter.   Similiarly, 
^^_  can  be  fixed  while  Cu   varies  from  zero  to  infinity  to  map  lines  of 
constant  z-plane  angle  for  the  complex  roots  onto  the  parameter  plane. 
While  the  loci  of  this  constant  angle  on  the  parameter  plane  does  not  have 
any  direct  usefulness  comparable  to  the  loci  of  constant  magnitude  and  the 
determination  of  settling  time,  they  do  lead  to  a  very  important  secondary 
property.   With  loci  of  constant  C5    and  constant  ^'  drawn  on  the  para- 
meter plane,  every  point  on  the  z-plane  now  has  a  unique  point  defined  on 
the  parameter  plane  in  the  case  of  the  second  order  systems.  The  fact  that 
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for  higher  order  systems,  several  z-plane  points  may  map  onto  one  para- 
meter plane  point  will  be  discussed  later. 

For  the  second  order  case,  this  one-to-one  correspondence  between 
points  on  the  z-plane  and  points  on  the  parameter  plane  is  useful  in  con- 
formal  mapping  of  various  loci  from  the  z-plane  to  the  parameter  plane. 
Among  the  possibilities  are  the  logarithmic  constant  damping  spirals  corres- 
ponding to  lines  of  constant  damping  in  the  s-plane.   Another  use  would  be 
the  mapping  of  lines  of  constant  overshoot,  [2]  and  [7],  from  the  z-plane 
to  the  parameter  plane.   This  z-plane  mapping  permits  the  designer  to 
select  feedback  compensation  which  simultaneously  satisfies  overshoot  and 
settling  time  criteria  for  the  second  order  system. 

2.3  Constant  Real  Root  Loci . 

In  order  to  completely  specify  all  possible  z-plane  root  locations  on 
the  parameter  plane,  real  root  configurations  must  be  allowed  for  as  well 
as  the  complex  root  configurations.   As  indicated  by  equation  (12),  the 
loci  of  constant  real  roots  are  straight  lines  on  the  parameter  plane. 
For  the  second  order  system,  the  closed  loop  poles  are  either  both  real  or 
both  complex.   Therefpre,  the  zone  of  real  root  description  and  the  zone  of 
complex  root  description  cannot  overlap,  except  at  a  critical  point  where 
the  real  root  is  about  to  emerge  from  the  real  axis. 

This  zone  separation  characteristic  greatly  simplifies  the  sketching 
of  the  parameter  plane  for  the  second  order  system.   Fortunately,  only  one 
set  of  calculations  is  necessary.   The  constant  0uy    curves  must  be  cal- 
culated.  However,  since  these  loci  of  constant  magnitude  are  straight 
lines,  only  a  minimum  amount  of  calculations  is  involved.   In  plotting 
the  constant   £°   curves,  one  has  only  to  connect  the  coordinates  of 
equal  ^   on  the   _\  curves.   Then,  to  plot  the  lines  of  constant  r?*~  , 
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since  from  equation  (2)  we  can  write, 

£rj=  <^X  5*  (24) 

one  has  only  to  draw  straight  lines  tangent  to  the  |  N  |  =  1.0  line  at 
the  t-3  intersections.   For  example,  a  line  drawn  tangent  to  the  N  =  -1.0 
line  at  the  intersection  of  the    ^^  =0.5  line  describes  the  real  root 
locations  (J^L    =   -0.5.   This  sketching  technique  will  not  necessarily  be 
extendable  to  the  higher  order  systems  unfortunately;  however,  it  does 
increase  the  simplicity  of  the  second  order  parameter  plane. 

Figure  (5)  is  the  parameter  plane  for  the  pure  second  order  system 
of  figure  2  which  now  has  some  of  the  loci  of  constant  angle  (solid  lines) 
and  some  of  the  loci  of  constant  real  roots  (broken  lines)  superimposed 
upon  it.  Note  that  this  parameter  plane  configuration  also  accurately 
describes  the  system  stability  in  terms  of  the  two  parameters. 


2.4  Pole  Effect. 

In  order  to  study  the  pole  effect  of  the  second  order  systems  on  the 
parameter  plane,  the  system  shown  in  figure  (6)  was  studied.  This  system 
can  also  be  represented  by  the  vector-matrix  equation  (14) ,  where  in  this 
case 
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Figure  6 
General   Second  Order 
System  Diagrams 
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A   = 
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The  characteristic  equation  is 

2 


-  A  z  +  AQ  ■  0 


(25) 


where 


and 


Aj  =  Ka^T/p-l/p+e'^/p2)  +  Ka2(l/p-e"pT/p) -1-e' 


PT 


AQ  =  Ka1(l/p2-Te"pT/p-e"pT/p2)  +  Ka2(e"pT/p-l/p)+e'pT. 
One  important  difference  between  the  pure  second  order  system  and  the 
system  with  damping  becomes  evident  by  inspecting  the  characteristic  equa- 
tion above.   In  the  pure  second  order  case,  by  judicious  definition  of  the 
system  parameters,  equation  (23),  it  was  possible  to  produce  a  parameter 
plane  as  a  function  (indirectly)  of  all  physical  parameters,  K,  T,  a.,  and 
a_.  By  moving  one  of  the  poles  into  the  left-half  of  the  s-plane,  it  is 
no  longer  possible  to  obtain  parameters  involving  T.   However,  the  system 
gain,  K,  remains  a  common  factor  in  the  characteristic  equation  coeffici- 
ents. Therefore,  for  the  second  order  system  with  damping,  we  can  define 
the  parameters  as 


C<   =  Ka 


1 


(26) 


A 


Ka, 


For  further  study  of  the  pole  effects,  let  p  ■  T 
tion  (9), 


1.0,  then  from  equa- 


d  »  0.368 

b  -  0.264 

c  -  -0.632 

dv   «  -J .368 

b]  -  0.368 

cx   =  0.632 

d2   -  1.0 

b2  =  o 

c2-0 
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and  from  equations  (10) , 

B1   =  -0.264U-1<  <^)     B2  =  0.368  ^U^  <^  ) 

Cj  =  0.632U_1(  J  )      C2  =  0.632  ^U^  f  ) 

Dj  -  -0.368U_1(  ^  )  -  ^  V  ^) 

D2=  -1.368  0)^(5;)+  &£u?<  5^). 
Figure  (7)  is  the  parameter  plane  for  this  system,  with  loci  plotted 
for  comparison  with  figure  (5).   Note  that  the  shifting  of  one  pole  from 
the  origin  of  the  s-plane  into  the  left-half  plane  has  produced  a  counter- 
clockwise rotation  and  a  compression  of  the  upper  semi-circle  contours. 
This  suggests  that  the  further  into  the  left  half  plane  of  the  real  pole 
(closer  to  the  origin  of  the  z-plane)  the  more  sensitive  the  system  re- 
sponse is  to  changes  in  a?.   Additionally,  from  figure  (7),  the  extremities 
of  the  complex  root  region  on  the  parameter  plane  have  been  extended,  to 
the  point  where  negative  values  of  beta  are  now  within  the  stable  region. 
From  figure  (7),  minimum  prototype  response  to  a  step  input  should  be  ex- 
pected for 

aC    =1.58 
fi     =  1.24 
Consequently,  there  are  two  ways  to  achieve  minimum  prototype  compensa- 
tion.  First,  if  K  is  fixed  to  some  value  other  than  1.58,  for  example 
1.0,  then  non-unity  position  feedback  must  be  used.  With  K  =  1.0,  the 
magnitudes  of  the  state  variables  at  consecutive  sampling  intervals  in 
response  to  a  unit  step  input,  calculated  using  equation  (14),  are  shown 
by  table  3.   Note  that  this  response  is  identical  to  that  shown  on  figure 
(3) ,  except  for  a  constant  scale  factor,  which  is  naturally  equal  to 
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TABLE  3 


nT 

0 
1 
2 
3 
4 


xt(t) 


0 

0.368 

0.632 

0.632 

0.632 


x2(t) 


0 

0.632 

0 

0 

0 


For  the  particular  case  where  the  conventional  error  definition  is 
required,  to  obtain  the  minimum  prototype  response,  K  must  be  restricted 
to  the  value  of  1.58.   Hence,  the  corresponding  value  for  a_  is  1.24/K  or 
0.785.   The  calculated  response  using  these  parameter  settings  is  shown 
by  table  4,  assuming  a  unit  step  input  again. 

TABLE  4 
nT 


0 
1 
2 
3 
4 


xt(t) 


0 

0.582 

1.000 

1.000 

1.000 


x2(t) 


0 

1.0 

0 

0 

0 


2.5  Zero  Effect. 

To  study  the  effect  of  a  real  zero  on  the  second  order  parameter 
plane,  the  system  shown  in  figure  (8)  was  analyzed.  For  this  system, 
the  state  transition  matrices  are: 


$  = 


1  -  Ka1(T  +  ZT  /2) 


Ka  (1  +  ZT) 


T  -  Ka2(T  +  ZT  /2) 


1  -  Ka2(l  +  ZT) 


and 


A    = 


K(T  +  ZT  /2) 
K(l  +  ZT) 
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Figure  8 
Second  Order  System 
with  Zero, 
Diagrams 
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(b)  Signal  Plow-Graph 
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Again,  the  characteristic  equation  is 


2 

z  +  Axz  +  A  =  0 


where 


AQ  =  Ka  (ZT2  -  2T)  -  Ka  (2ZT  +  2)  +  1 

2 
and  A  =  Ka^ZT  +  2T)  +  Ka2(2ZT  +  2)  -  2 

Again,  a  choice  of  parameters  for  the  parameter  plane  can  be  made  as  was 
performed  in  section  2.4.   For  further  study,  Z  and  T  must  be  fixed  in 
magnitude.   For  this  example  then,  let  Z  ■  T  =  1.0.   Following  the  now 
familiar  calculations,  the  parameter  plane  of  figure  (9)  results.   Similiar 
to  the  results  shown  for  the  pole  effect,  the  parameter  plane  loci  have 
been  rotated  counter-clockwise  from  the  pure  second  order  configuration  of 
figure  (5).   Apparently  this  counter-clockwise  rotation  can  be  associated 
with  the  relative  stability  of  the  system,  as  both  the  system  with  the  zero 
and  the  system  with  the  left-half  s-plane  pole  are  inherently  more  stable 
than  the  pure  second  order  servo.   However,  no  measure  for  quantitatively 
describing  this  relative  stability  has  been  discovered. 

From  figure  (9),  minimum  prototype,  ripple  free  response  to  a  step  in- 
put is  predicted  for 

c<   -  0.5 
/&     -  0.125 
Again  there  are  two  choices  for  the  physical  parameters  in  achieving 
this  compensation.   If  K  is  variable,  unity  position  feedback  may  be  used. 
If  K  is  fixed,  then  non-unity  position  feedback  is  required.   Responses 
similiar  to  that  achieved  in  the  last  section  result  for  each  of  these 
options. 
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CHAPTER  111 
HIGHER  ORDER  SYSTEMS 

3.1  Parameter  Limitations. 

It  can  be  readily  demonstrated  that  in  order  to  achieve  complete  free- 
dom in  locating  the  closed  loop  roots  of  an  n-th  order  characteristic 
equation,  the  control  system  error  signal  must  contain  feedback  intelli- 
gence on  each  of  the  n  dynamically  independent  "states"  of  the  system. 
Consequently,  the  system  analog  or  mathematical  model  consists  of  n  para- 
meters and  the  parameter  space  becomes  an  n-space.   Unfortunately,  the 
existing  Mitrovic,  or  algebraic,  methods  have  their  basis  in  the  solution 
of  two  simultaneous  equations,  equations  (5a)  and  (5b).   This  restricts 
the  mathematical  treatment  of  any  n-th  order  system  to  only  two  degrees  of 
freedom.  While  this  reduction  in  the  degrees  of  freedom  restricts  the 
ability  to  position  the  closed  loop  poles  in  the  z-plane,  a  large  class 
of  control  systems  can  be  reduced  to  the  two  parameter  problem. 

Typically,  the  gain  setting,  K,  in  figure  (10)  is  dictated  by  the 
requirements  for  steady  state  error.   A  second  degree  of  freedom  is  normal- 
ly eliminated  by  the  adoption  of  the  conventional  error  definition  which 
demands  unity  feedback.   Finally,  in  most  physical  systems ,  noise  free 
derivative  signals  higher  than  the  second  derivative  are  unobtainable. 
Thus,  the  only  degrees  of  freedom  for  feedback  compensation  typically  avail- 
able are  the  magnitudes  of  the  first  and  second  derivative  feedback. 

With  this  restriction  on  the  positioning  of  the  closed  loop  roots  in 
mind  then,  the  feedback  compensation  of  systems  with  order  greater  than  two 
can  be  designed  on  the  parameter  plane  using  techniques  already  discussed. 
Two  additional  features  of  the  parameter  plane  must  now  be  considered. 
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Figure  10 
Typical  Block  Diagram  for 
n-th  Order  Sampled- Data  System 
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Unfortunately,  at  the  present  time,  these  two  features  tend  to  reduce 
the  effectiveness  of  the  Mitrovic  Method*   However,  it  is  believed  that 
this  reduced  effectiveness  is  still  an  improvement  over  the  existing  one 
parameter  design  techniques. 

3.2  Maximum  Modulus  Uncertainity . 

The  parametric  equations  describing  the  parameter  plane  loci  were  de- 
rived from  the  system  characteristic  equation.   Any  specified  location  (M- 
point)  on  the  resulting  parameter  plane  can  be  interpreted  as  the  values 
of  the  two  parameters  which  will  put  a  closed  loop  root  at  the  specified 
position  of  the  z-plane.   Therefore,  for  a  typical  fourth  order  system 
for  example,  the  closed  loop  root  configurations  shown  by  figures(lla)  and 
(lib)  each  correspond  to  a  unique  point  on  the  parameter  plane,  each  lying 
on  the  (k).,   =  1.0  loci.   In  order  to  distinguish  between  the  different  con- 
figurations on  the  parameter  plane,  it  is  also  necessary  to  plot  the  CJ '-   0.5 
and  Ck),=   1.5  loci.   The  intersection  of  two  of  the  loci  of  constant  ^3^ 
then  adequately  describes  the  location  of  each  pair  of  complex  roots.   If 
figure  (12)  represents  the  parameter  plane  for  the  fourth  order  system  of 
figure  (11),  then  the  root  location  depicted  by  figure  (11a)  corresponds 
to  point  A,  while  the  configuration  of  figure  (lib)  corresponds  to  point 
B  on  the  parameter  plane. 

Thus,  for  an  n-th  order  system,  the  modulus  of  each  real  root  or  each 
pair  of  complex  roots  is  completely  determined  on  the  parameter  plane. 
The  accuracy  of  that  determination,  however,,  is  directly  related  to  the 
number  of  loci  of  constant    ^   drawn.   In  the  second  order  case,  it  was 
fairly  easy  to  extrapolate  between  adjacent  loci.   As  the  order  of  the 
characteristic  equation  increases,  the  difficulty  of  making  this  extra- 
polation increases.   It  should  be  noted  that  this  problem  of  interpreting 
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Figure  11 
Typical  z-Plane  Closed  Loop  Root 
Configuration  for 
Fourth  Order  System 
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intersections  of  constant  loci  also  exists  in  the  second  order  case  when 
the  pl^nt  gain  is  insufficient  to  produce  complex  roots.   Uncertainity 
arises  as  the  order  of  the  system  increases.   For  example,  given  a  fifth 
order  system,  and  a  corresponding  parameter  plane  M-point  defined  by  the 
intersection  of  two  constant  ^L  loci  and  one  real  root  loci.   Obviously, 
this  M  point  must  determine  two  pair  of  complex  roots  and  one  real  root. 
However,  with  existing  techniques,  it  is  impossible  to  associate  any 
specific  Oj     with  a  specific  complex  root. 

3.3  Discontinuities  on  the  Parameter  Plane . 

A  second  feature  which  may  exist  in  higher  order  systems  involves  dis- 
continuities on  the  constant  CJ     and   C   loci,   Equation  (11)  describes 
the  parametric  equations  for  alpha  and  beta  as  a  ratio  of  two  polynomials. 
Consequently,  whenever 

31C2  '  B2C1  =  °'  (26) 

a  discontinuity  exists  on  the  loci  of  constant  Cj    or  C"   .At  the  present 

time,  the  physical,  significance  of  this  discontinuity  is  not  understood, 

in  fact  one  may  not  even  exist.   Its  presence,  however,  does  complicate 

the  sketching  and  interpretation  of  the  parameter  plane. 

3.4  Dominance  and  Settling  Time . 

In  the  case  of  the  linear,  continuous  system,  the  response  of  a 
higher  order  system  can  often  be  discussed  in  terms  of  a  dominant  pair  of 
roots  in  the  s-plane.   This  is  a  pair  of  complex  conjugate  roots  whose 
residue  and/or  transient  settling  time  leads  to  by  far  the  greatest 
portion  of  the  observed  physical  response.   In  the  s-plane,  these  roots 
have  an  argument  and  a  real  part  considerably  less  in  magnitude  than  any 
other  pair  of  complex,  or  real,  roots.   If  such  a  pair  of  roots  exist,  the 
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responnc  can  be  quite  accurately  predicted  in  terms  of  a  second  order 

approximation  and  the  settling  time  predicted  by  equation  (17), 

Similiarly,  for  the  sampled-data  system,  the  dominance  of  the  system 

response  can  be  discussed  in  terms  of  the  pair  of  complex  roots,  or  real 

root,  whose  modulus  or  argument  is  considerably  greater  than  that  of  any 

other  root,  real  or  complex.   If  such  a  root  exists,  the  settling  time 

can  be  closely  approximated  by  a  slight  modification  to  equation  (22)  such 

that 

t  =  nT  -  4T/ln  (    )        .  (27) 

s  max 

When  the  dominance  cannot  be  clearly  determined  by  the  difference  in  the 
various  arguments,  equation  (27)  becomes  less  and  less  accurate  in  pre- 
dicting the  system  settline  time. 

Alternatively,  the  minimization  of  the  settling  time  can  be  achieved 
qualitatively  by  minimizing  the  modulus  of  each  root  or  pair  of  complex 
roots.   That  is  to  say,  compensation  must  be  chosen  which  will  place  all 
the  z-plane  roots  within  a  circle  of  minimum  radius.   If  this  circle  of 
minimum  radius  happens  to  be  the  origin  of  the  z-plane,  then  minimum  pro- 
totype response  is  guaranteed.   If,  however,  all  the  z-plane  roots  cannot 
be  placed  at  the  origin  due  to  the  limitations  on  the  degrees  of  freedom 
available,  minimum  settling  time  will  be  somewhat  greater  than  minimum 
prototype,  and  will  not  be  ripple  free.   On  the  parameter  plane,  this 
minimization  of  the  settling  time  involves  the  examination  of  the  inter- 
sections of  the  constant  V   loci  for  each  permissible  M-point,  and 
selecting  that  one  which  does  place  all  the  roots  within  the  circle  of 
minimum  possible  radius.   Observations  of  the  resulting  arguments  will 
enable  the  designer  to  determine  if  equation  (27)  is  then  valid  or  not. 
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3 .5  Third  Order  Example . 

Consider  the  system  shown  on  figure  (13).   This  system  can  be  re- 
presented by  the  signal  flow  graph  of  figure  (14),  and  the  vector-matrix 
equation 

X  [(k  +  1)t]  =  §    2L<kT)  +  ^  e<kT) 


where 


J, 


-T        -?T  -T        -2T 

(3/2  -  2e   +  l/2e  l)      (1/2  -  e   +  l/2e  L) 


-T    -2T 
(2e   -  e  Zi) 

-T      -2T 
(-2e  l   +  2e  *L) 


(e   -  e   ) 

-T     -2T 
(-e  l  +   2e   ) 


A  = 


T/2  -  3/4  +  e"T  -  l/4(e"2T) 

-T         -2T 
1/2  -  e  l   +  l/2(e  L) 

-T    -2T 
e   -  e 


and,  e(kT)  =  R(kT) -X1(kT) -a1X2(kT) -a2X3(kT) 

The  characteristic  equation,  evaluated  at  T  =  1  second  is 

3      2 

Z  +  A2Z  +  A  Z  +  AQ  =  0 


with 


A2  =  0.1998a  +  0.2325a2  -  1.4191 
A  =  -0.1263a  -  0.4650a  +  0.7238 


and  A  =  -0.0734a  +  0.2325a  -  0.0310 

Observe  that  it  is  no  longer  possible,  as  was  done  for  the  second  order 
systems,  to  include  the  plant  gain,  K,  in  the  definition  of  the  system 
parameters.   This  combination  of  the  gain,  K,  with  the  feedback  coeffici- 
ents is  only  possible  so  long  as  each  coefficient  in  the  numerator  of  the 
product  G(s)H(s)  is  some  product  containing  a  variable  feedback  gain.   In 
the  second  order  case, 

G(s)H(s)  =  K(ax  +  a2s) 


D(s) 

The  characteristic   equation  can  be   formed  as 

1  +  z-transform  of      Fg(s)H(s)J 
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Thus,  the  only  place  that  the  feedback  coefficients  occur  in  the  character- 
istic equation,  they  are  each  multiplied  by  the  gain,  K.   However,  for 
this  third  order  example,  the  produce  G(s)H(s)  has  the  form, 

G(s)H(s)  =  K(  1  +  a^  +  a?s2) 


D  (s) 
In  the  characteristic  equation  then,  K  is  a  multiplier  of  terms  other  than 
those  involving  the  feedback  coefficients,  a.,  and  a„.   For  the  example,  we 
are  restricted  to  defining 


fi   s* 


A 

a 

2 


Following  the  usual  manipulations,  the  third  order  parameter  plane,  figure 
(15),  was  achieved.   No  lines  of  constant   V  were  included  on  the  para- 
meter plane,  only  to  avoid  unnecessary  cluttering.   Observe  that  the  separa- 
tion of  complex  root  regions  and  real  root  regions  no  longer  exists;  which 
should  be  expected,  as  the  root  locus  always  has  one  real  root.   Also  note 
that  the  region  of  stability  is  that  region  of  the  parameter  plane  which 
is  enclosed  by  the  contours  Cj  =  1.0,   0^  =  -1.0,  <J^   =  1.0.   From  figure 
(15),  minimum  settling  time  using  feedback  compensation  is  predicted  when 

ax  =  1.61 

a2  =  0.54 
The  circle  of  minimum  radius  enclosing  the  roots  is  ^L-   0.3. 
Table  5  lists  the  calculated  magnitudes  of  the  state  variables  at  consecu- 
tive sampling  intervals  in  response  to  a  unit  step  input  using  these  values 
for  the  feedback  coefficients.   From  Equation  (27)  then,  minimum  settling 
time  is  calculated  as 

t  =  3T  -  4T/ln  0.3 
s 

-  6.3T 
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At  time  7T,  table  5  shows  that  the  output  has  settled  to  within  57„  of  the 
final  value.   Figure  16  shows  the  state  variable  time  response  for  this 
minimum  settling  time  compensation.   These  results  verify  the  effective- 
ness of  equation  (27)  in  predicting  the  servo  settling  time  to  a  good 
engineering  approximation,  even  with  the  dominance  of  the  roots  at  maxi- 
mum modulus,  as  discussed  in  section  3.4,  somewhat  uncertain. 

TABLE  5 


nT 

xt(t) 

0 

0 

1 

0.08405 

2 

0.33610 

3 

0.57736 

4 

0.74605 

5 

0.85111 

6 

0.91357 

7 

0.95002 

8 

0  97114 

9 

0.98334 

10 

0.99039 

11 

0.99445 

12 

0.99680 

13 

0.99815 

14 

0.99893 

15 

0.99938 

16 

0.99965 

17 

0.99980 

18 

0.99988 

19 

0.99993 

20 

0.99999 

x2(t) 


0 

0.19979 

0.26768 

0.20644 

0.13340 

0.08048 

0.04724 

0.02743 

0.01586 

0.00916 

0.00529 

0.00305 

0.00176 

0.00102 

0.00059 

0.00034 

0.00020 

0.00011 

0.00007 

0.00004 

0.00001 


x3(t) 


0 

0.23254 
-0.00653 
-0.06887 
-0.05968 
-0.03964 
-0.02411 
-0.01419 
-0.00824 
-0.00477 
-0.00275 
-0.00159 
-0.00092 
-0.00053 
-0.00031 
-0.00018 
-0.00010 
-0.00006 
-0.00003 
-0.00002 
-0.00001 


3.6  Algebraic  Design  Methods. 

One  simple  and  quite  obvious  method  for  designing  feedback  compensa- 
tion for  minimum  prototype  response  should  not  be  overlooked.   In  determin- 
ing the  system  characteristic  equation  as  a  function  of  the  feedback  vari- 
ables, which  must  be  done  to  apply  the  Mitrovic  Method,  an  alternate  method 
arises  for  determining  minimum  prototype  compensation.   To  obtain  minimum 
prototype  response  for  a  given  system,  the  function  of  any  compensation 
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Figure  16 
State  Variable  Time  Response  for  Parameter 
Plane  Compensated  Third  Order  System 
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must  be  to  drive  all  the  closed  loop  roots  to  the  origin  of  the  z-plane. 

The  characteristic  equation  of  an  n-th  order  system,  such  as  shown  on 

figure  (10),  is 

zn  +  A  ,zn4  +  +  An  =  0  (28) 

n-1  0 


where 


A.  =  B_  .  a„  +  B.  .  a.    + +  B  .  a  +  B 

x    0i  0    li  1  ni  n    n 


and  a.  =  feedback  coefficient  of  the  i-th  derivative 

J 

feedback.   Therefore,  to  algebraically  design  the  minimum  prototype  compen- 
sation, thereby  positioning  each  closed  loop  root  at  the  origin  of  the  z- 
plane,  one  has  only  to  simultaneously  solve  the  n  non-homogeneous  equa- 
tions 

Vi  ■  ° 

A„-2  "  ° 

The  solution  to  these  simultaneous  equations  determines  the  values  of  the 
n  feedback  parameters  required  to  produce  minimum  prototype  response. 
Additionally,  as  was  pointed  out  in  section  3.5,  since  the  number  of  feed- 
back variables  now  equals  the  degree  of  the  characteristic  equation,  we 
can  rewrite  the  A.  in  equation  (28)  as 

A.  =  Bn.Ka„  +  B,  .Ka,  + +  B  .Ka  +  B       , 

i    Oi  0    li  1  ni  n    n 

thereby  providing  the  designer  with  yet  another  potential  variable. 

To  demonstrate  the  use  of  this  algebraic  design  technique,  consider 
again  the  third  order  example  of  figure  (13) .   Let  the  open  loop  gain  be 
variable,  K,  instead  of  1.0  as  shown.   Furthermore,  let  the  coefficient  of 
the  position  feedback  path  be  a   .      Then,  for  minimum  prototype  response, 
the  characteristic  equation  becomes 
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3      2 
z  +  A  z  +  A.z  +  A  -  0 

where     A  =  0.2325Ka2  +  0.1998Ka  +  0.0840Ka  -  1.503  =  0 

A  =  -0.4650Ka  -  0.1263Ka  +  0.1709Ka  +  0.5529  =  0 

A„  =  0.2325Ka„  -  0.0734Ka,  +  0.0188Kart  -  0.0498  -  0 
0  2  10 

The  simultaneous  solution  of  these  equations  then,  is 

KaQ  =  3.652 

Ka  -  4.448 

Ka  =  1.323 
In  order  to  permit  the  adoption  of  the  conventional  error  definition, 
then,  K  =  3.652 

aQ  =  1.0 

a  =  1.218 

a2  =  0.363 
Using  these  values  of  the  system  parameters,  the  response  of  the  system 
of  figure  (10)  to  a  unit  step  input  is  as  shown  by  table  6.   The  time 
response  of  the  state  variables  is  shown  on  figure  (17). 

TABLE  6 


nT 

xt(t) 

x2(t) 

x3(t) 

0 

0,000 

0.000 

0.000 

1 

0.307 

0.730 

0.849 

2 

0.929 

0.286 

-0.849 

3 

1,000 

0.000 

0.000 

4 

1.000 

0.000 

0.000 
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Figure   17 
State  Variable  Time  Response  for  Algebraic 
Compensated  Third  Order  System 
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CHAPTER  IV 
DIGITAL  COMPUTATION  IN  MULTI- ORDERED  PARAMETER  SPACE 

While  the  feedback  compensation  of  many  n-th  order  sampled-data  con- 
trol systems  can  be  reduced  to  the  two  parameter  problem  presented  in 
Chapter  III,  circumstances  can,  and  do,  arise  where  derivative  signals 
higher  than  the  second  derivative  are  available.   Since  each  additional 
degree  of  freedom  available  to  the  designer  enlarges  the  region  in  the  z- 
plane  in  which  the  closed  loop  roots  may  be  positioned,  common-sense  alone 
dictates  their  usage  where  available.   As  pointed  out  in  preceding  chapters, 
Mitrovic's  Method  is  primarily  two  parameter  in  nature.   Additional  design 
techniques  are  then  desirable  whereby  a  parameter  space  higher  than  a  2- 
space  can  be  achieved. 

The  digital  computer  provides  one  method  of  increasing  the  order  of 
the  parameter  space.   Any  design  process  or  technique  is  simply  a  method 
for  locating  the  roots  of  a  polynomial,  the  system  characteristic  equation, 
as  a  function  of  one  or  more  parameters.   These  methods  were  developed  pri- 
marily to  surmount  the  difficulty  of  factoring  a  given  polynomial.  With 
the  digital  computer  however,  this  polynomial  factoring  can  be  repeatedly 
performed,  always  with  greater  accuracy,  and  sometimes  in  less  total  time 
than  analysis  performed  by  the  graphical  and  semi-graphical  techniques.   In 
this  chapter,  one  approach  to  multi-parameter  study  using  the  digital  com- 
puter will  be  presented.   The  emphasis  on  the  application  of  this  method 
will  be  directed  toward  feedback  compensation;  however,  the  extension  to 
cascade  compensation  or  other  parameter  variations  does  exist. 

Recent  emphasis  on  the  state  space  method  of  analysis,  [2],  as  applied 
to  sampled-data  systems  affords  an  excellent,  if  not  necessary,  point  of 
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departure  for  digital  computer  analysis.  The  general  form  of  the  state 
space  vector-matrix  equation  describing  a  sampled-data  feedback  control 
system  was  given  by  equation  (14) , 

_X  [(k  +  1)T]  =  J  X(kT)  +  A    R(kT) 
It  can  easily  be  demonstrated,  see  appendix  1,  that  the  characteristic 
equation  of  the  state  transition  matrix,  jj   is  indeed  the  characteristic 
equation  of  the  closed-loop  transfer  function.   Since  numerous  numerical 
methods  presently  exist  for  determining  the  eigenvalues  of  a  given  matrix, 
the  determination  of  the  closed  loop  pole  locations  by  evaluating  the 
eigenvalues  of  the  jf     matrix  is  readily  performed  on  the  digital  computer. 
To  determine  these  root  locations  as  a  function  of  one  or  more  parameters, 
one  has  only  to  repeatedly  increment  the  parameters  and  determine  the  eigen- 
values of  the  resulting  M    matrix.   In  this  fashion,  an  n-coordinate 
space  can  readily  be  created  by  computer  programming  techniques.   One 
method  is  to  use  the  FORTRAN  language,  and  program  DO  LOOP's  inside  other 
DO  LOOP's  until  sufficient  coordinates  are  defined  to  permit  the  designer 
to  sweep  through  a  desired  region  in  the  n-space,  determining  the  root 
locations  corresponding  to  each  point  in  the  n  dimensional  parameter  space. 

In  order  to  demonstrate  this  method  of  analysis  more  fully,  the 
following  specific  project  was  undertaken.   Given  the  perfectly  general 
second  order  sampled-data  system  shown  in  figure  (18),  the  values  of  the 
parameters  Ka  and  Ka  required  for  minimum  prototype,  minimum  settling  time 
response  to  a  step  input  are  to  be  determined  as  functions  of  the  sampling 
interval,  T.   The  resulting  computer  program  should  be  such  that  the  de- 
sired design  data  can  be  readily  generated  with  only  a  minimum  amount  of 
additional  input  data  supplied  to  the  data,  for  any  second  order  system 
with  real  poles. 
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The  system  of  figure  (18)  can  be  represented  by  the  open-loop  vector 
matrix  equation: 

X[(k  +  1)T]  =  FX(kT)  +  A  e(kT). 
The  error  signal  at  the  k-th  sampling  interval  is  represented  as; 

e(kT)  =  R(kT)  -  a^CkT)  -  a2X2(kT) 
In  the  chosen  example,  the  F  matrix  is  perfectly  general,  irregardless  of 
the  location  of  the  olant  poles  and  zeros,  and 


F  -  Q, 


-p.T     -p0T 
p  e  '  1   -p  e  v2 


0,(^1  +  °2>T  -  i)/o0 


e"PlT  -  e"^ 


-Ple-PlT  +  p2e-V 


where 


0 


l/(n2-P;L) 


°1  =  P1P2/(P!  +  P2> 
Modifications  to  the  F  matrix  are  needed,  however,  if  a  pure  second  order 

system  is  used.   The  &    matrix  on  the  other  hand,  is  dependent  upon  the 

number  of  oure  integration  s  in  the  forward  path  of  the  plant  transfer 

function.   For  example,  for  p   =  p_  =  0, 


A  =K 


b2  +  bLT  +  bQT  II 


*>2d<T>   +  bx  +  bQT 
When  only  one  pure   integration  exists,    that   is   p.  ^:0,      p     =  0, 
(1/Ptl)2    [(bj  +  b0T)Pl    -   bQ  +  02e"PlT] 

b2   I(T)   +  VP1    -   V'PlT/pl 


A=   K 


where, 


°2      =      (b2Pl    "  VP1   +  b0 
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General  Second  Order 
System  With  Zeros 


Finally,   when   n     ^    0,    and   P9  ^    0, 


A^  K 


where 


VP1P2  +  °2°3e"PlT  +  W"'21 


b2^"(T)  +  G1Q1e"PlT  -G2Oie"P2T 


a     2 
Q3  -  l/(Pj_  -  P]P2y> 

°4  ^  X^P2  "  P1*V 


Gl  "  Pl  "  Plbl  +  bo 
G2  *  P2  "  p2bi  +  b0 


o5  .  (b2p2  -  bl)p2  +  bQ 


d(T)  =  unit  impulse. 
The  program  which  was  developed,  see  appendix  2,  to  satisfy  the 
design  goal  requires  four  inputs  in  the  form  of  data.   In  the  order  in 
which  they  are  required,  they  are: 

a.  the  pole  locations  of  the  plant  in  the  s-plane.   The  nega- 
tive value  of  this  s-plane  location  must  be  entered  on  the  data  card,  and 
the  poles  should  be  listed  in  decreasing  magnitude  order. 

b.  the  plant  s-plane  numerator  coefficients  listed  in  ascend- 
ing order. 

c.  the  minimum  and  maximum  values  of  T,  the  sampling  interval, 
and  the  incremental  value  of  the  sampling  interval  to  be  used  in  the  in- 
vestigation.  T  minimum  must  be  greater  than  zero. 

d.  miscellaneous  title  information  and  graph  scales  required 
for  graph  output  in  accordance  with  an  existing  graph  subroutine  (DRAW) . 

The  program  outputs  are: 

a.  a  tabulated  listing  of  the  minimum  settling  time  values  of 
Ka-  and  Ka  determined  for  each  value  of  T  used  in  the  investigation. 

b.  a  graph  plot  of  the  minimum  settling  time  value  of  Ka..  vs. 
T. 

c.  a  graph  plot  of  the  minimum  settling  time  value  of  Ka-  vs. 


48 


Let  us  now  consider  the  geometric  significance  of  this  program.   It 
has  been  demonstrated  that  for  a  given  value  of  the  sampling  interval  T, 
the  origin  of  the  z-plane  maps  into  a  single  point  on  the  parameter  plane 
for  a  second  order  system.   However,  if  a  three  dimensional  right-handed 
coordinate  set  is  created,  whose  axes  are  Ka.  ,  Ka_,  and  T,  the  minimum 
settling  time  point  corresponding  to  6»X  =  0  will  map  into  a  curve  in  the 
3-space.   The  plots  of  Ka  and  Ka»  required  for  minimum  settling  time  re- 
sponse vs.  T  can  be  thought  of  as  the  projection  of  this  curve  in  3-space 
onto  the  Ka  -T  and  Ka  -T  planes.   Similiarly,  any  n-th  order  system,  with 
m  degrees  of  freedom,  where  m  ^  n,  can  be  investigated  in  a  computer 
oriented  m-space  and  displayed  to  the  designer  by  a  set  of  m  projections 
on  the  a.  vs.  T,  (i  =  1,  2,...,m),  planes.   If  m  =  n,  a  set  of  projections 
on  the  Ka.-T  planes,  (i  =  1,  2,...,n),  are  available. 

The  following  examples  will  serve  to  illustrate  the  class  of  feed- 
back compensation  problems,  the  solutions  of  which  can  be  greatly  simpli- 
fied by  using  this  particular  programming  technique. 

EKAMPLE  #1 

A  sampled-data  control  system  whose  plant  transfer  function  is 

G(s)  =  K(s  +  l)/s(s  +  2) 
is  to  be  used  with  a  dual  capacity  for  sampling  rates  of  lOcps  and  leps. 
Compensation  is  to  be  designed  to  provide  minimum  prototype  response  to 
a  step  input  for  each  of  the  possible  sampling  frequencies. 

Figures  (19)  and  (20)  are  the  minimum  settling  time  parameter  planes 
for  this  plant  transfer  function.   From  these  parameter  planes,  the  re- 
quired parameter  values  for  each  sampling  frequency  can  be  obtained  as: 

T  Ka  Ka. 

0.1  sec  10.0  0.949 

1.0  sec  1.15  0.538 
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Finally,  for  unity  position  feedback,  the  physical  plant  parameters 

must  be: 

T  Ka  Ka. 

—1  —2 

0.1  sec  10.0  0.095 

1.0  sec  1.15  0.467 

Figure  (21)  represents  the  block  diagram  of  the  compensated  system. 
When  the  three  switches  in  the  physical  plant  are  "ganged",  then  the  control- 
led switch  of  the  sampling  frequency  will  introduce  the  correct,  associated 
feedback  compensation  for  minimum  settling  time. 
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EXAMPLE  #2 

A  system  whose  plant  transfer  function  is 

G(s)  =  K(s  +  Z)/s(s  4-  2) 
is  to  be  controlled  by  a  sampled  error  signal.   The  plant  zero  location 
is  a  function  of  a  cascaded  compensator,  and  can  be  placed  anywhere  on 
the  real  axis  of  the  s-plane.   A  sampling  interval  of  0.5  seconds  is 
intended.   The  designer  is  to  set  the  gain,  K,  the  zero  location,  Z,  and 
design  feedback  compensation  for  minimum  settling  time  in  response  to  a 
step  input,  and  minimum  steady  state  error  for  a  ramp  input. 

Using  the  program  to  determine  the  minimum  settling  time  values  of 
Ka.  and  Ka  as  a  function  of  T,  a  family  of  curves,  with  parameter  Z, 
can  be  constructed  on  the  Ka..  vs.  T  and  Ka.  vs.  T  planes.   These  families 
of  curves  are  shown  by  figures  (22)  and  (23) .   Of  these  two  families  of 
curves,  the  minimum  settling  time  Ka  vs.  T  plane  is  by  far  the  most 
interesting.   The  minimum  settling  time  Ka.  vs.  T  loci  was  shown  by  fig- 
ure (19)  to  be  hyperbolic.   Furthermore,  section  3.3  outlined  the  presence 
of  discontinuities  on  parameter  planes  for  systems  with  more  than  two 
parameters.   These  discontinuities  again  appear  on  the  minimum  settling  time 
parameter  planes,  figures  (22)  and  (23),  and  are  best  demonstrated  by  fig- 
ure (23).   Inspection  of  figure  (22),  the  minimum  settling  time  Ka.  vs.  T 
plane,  also  reveals  these  discontinuities  on  the  loci  of  minimum  settling 
time  for  the  cases  of  Z  less  than  zero  (right-half  s-plane  zeros). 

It  has  been  demonstrated,  [2]  and  [7],  that  for  a  type  one  sampled- 
data  control  system,  the  steady  state  error  in  response  to  a  ramp  input  is 
inversely  proportional  to  the  plant  root  locus  gain;  similiar  to  the  "velo- 
city lag"  error  in  continuous  systems.   Consequently,  to  minimize  the  ramp 
error  and  at  the  same  time  design  for  minimum  settling  time  response  to  a 
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step  input,  the  designer  must  maximize  the  parameter  K  (or  Ka  when  using 
unity  position  feedback)  on  the  specified  line  of  constant  sampling  in- 
terval T  =  0.5.   Inspection  of  figure  (22)  shows  the  importance  played  in 
this  selection  of  the  plant  gain  by  the  discontinuity  phenomena.   From 
figure  (22),  the  maximum  Ka  attainable  for  T  =  0.5  is  Ka  =  12.7,  when 
Z  =-1.5.   For  Z  less  than  -1.5,  the  discontinuity  asymptote  lies  to  the 
left  of  T  =  0.5,  and  a  negative  value  of  Ka..  (positive  feedback)  is  re- 
quired for  minimum  step  input  settling  time.   For  Z  greater  than  -1.5, 
lower  values  of  Ka.  are  required  for  minimum  settling  time  for  a  step  in- 
put and  therefore  a  larger  ramp  error  is  experienced. 

Therefore,  the  design  criteria  of  this  example  is  achieved  for: 

K  =  12.7 

ax   =   1.0 

a2  =  0.656 
Z  =  -1.5 
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CHAPTER  V 
CONCLUSIONS 

The  application  of  Mitrovic's  Method  to  feedback  compensation  of 
sampled-data  control  systems  has  been  shown  to  be  quite  useful  in  the 
minimization  of  system  settling  time.   The  parameter  plane  for  the  second 
order  system  presents  an  easily  obtainable  picture  of  the  system  stability 
and  settling  time  as  a  function  of  the  system  parameters.   Using  the 
methods  described  herein,  the  second  order  parameter  plane  development 
results  in  a  two  parameter  analysis  method  requiring  about  the  same  total 
effort  to  construct  as  conventional  one  parameter  methods. 

In  the  case  of  higher  order  systems,  the  effectiveness  of  the  Mitrovic 
Method  as  a  design  tool  is  somewhat  reduced,   Foremost9  the  plotting  of  the 
various  loci  requires  at  least  a  desk  calculators,  preferably  a  digital 
computer.   If  the  number  of  system  parameters  can  be  reduced  to  only  two, 
then  the  parameter  plane  is  still  a  good  tool  for  minimizing  settling  time. 
Unfortunately,  the  determination  of  the  minimum  settling  time  compensation 
from  the  higher  order  parameter  plane  requires  a  more  tedious  study  of  the 
plane  than  for  the  second  order  case.   The  intersections  of  all  the  loci 
of  constant  ^O^  and  constant  ($1     must  be  examined  in  detail  in  search  of 
the  circle  of  minimum  radius  in  the  z-plane  encircling  all  the  closed  loop 
roots. 

For  higher  order  systems  with  more  than  two  degrees  of  freedom, 
analysis  with  the  digital  computer  is  most  effective,  indeed  necessary. 
The  vector-matrix  state  space  equations  representing  the  given  system 
were  shown  to  lead  directly  to  the  determination  of  the  closed  loop  root 
positions  as  a  function  of  up  to  n  parameters,  n  being  the  order  of  the 
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system.   The  strength  of  this  method  is  the  ability  to  reduce  the  design 
criteria  to  a  function  of  almost  any  physical  parameter  or  condition. 
However,  since  a  large  fraction  of  the  feedback  compensation  problems  are 
not  more  than  two  parameter  in  nature,  the  Mitrovic  Method  can  be  an 
exceedingly  important  tool  to  the  design  engineer. 
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APPENDIX  I 
STABILITY  ANALYSIS  OF  VECTOR -MATRIX  DIFFERENCE  EQUATIONS 
The  general  form  of  the  state  space  vector-matrix  equation  describ- 
ing a  linear  system  with  a  sampled  input  is, 

X[(k  +  1)T]  =  FX(kT)  +  A   u(kT).  (1-1) 

In  this  form,  the  states  of  the  system,  X,  are  defined  as  the  control 
system  output  and  its  first  n-1  derivatives,  where  n  is  the  order  of  the 
linear  system.   For  a  feedback  control  system,  the  control  vector,  u(t), 
is  the  difference  between  the  command  input,  R(t),  and  the  sum  of  the 
feedback  quantities,  ie; 

u(t)  =  R(t)  -  AC  x(t),  (I_2) 

where  A  is  the  transpose  of  A,  the  coefficient  vector  controlling  the 
magnitude  of  feedback  for  each  state  variable.  In  discrete  time  form, 
equation  (12)  is  written  as 

u(kT)  =  R(kT)  -  Afc  X(kT)  .         (1-3) 

By  substituting  equation  (1-3)  into  equation  (1-1) ,  the  familiar  vector- 
matrix  difference  equation  for  a  sampled  data  control  system  is  obtained 
as 

X[(k  +  1)T]  =  FX(kT)-^  Afc  x(kT)  4-  A  R(kT)  ,     (1-4) 

or  by  imbeding  the  feedback  quantities  directly  into  the  F  matrix, 

X[(k  +  1)T]  =  J  X(kT)  +  A  R(kT)       .         (1-5) 
With  the  feedback  thus  imbeded  in  the  transition  matrix,  jf  ,  the 
stability  analysis  of  a  given  system  can  be  confined  to  an  analysis  of  the 
partial,  state  equation, 

X[(k  +  1)T]  =  _?X(kT)  ,         (1-6) 

which  is  an  initial  value  problem.   If  the  state  trajectory  returns  to  a 
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given  eouilihrium,  or  singular  point,  in  the  state  space  for  any  initial 
state  vector,  it  will  also  move  to  a  predictable  translation  of  this  singu- 
lar point  and  achieve  equilibrium  for  the  input  R(t)  which  can  be  described 
by  the  initial  state  vector.   The  dynamic  stability  of  the  system  is  thus 
ensured. 

Hahn,  [1],  studied  the  stability  of  equation  (1-6)  using  Liapunov's 
Direct  Method,  and  concluded  that  the  stability  of  equation  (I«6)  was 
ensured  when  the  logarithm  of  the  eigenvalues  of  the  <£  matrix  were  all 
less  than  zero  in  magnitude,  that  is 

In  Xi    <  °  •         d"7) 

This  is  equivalent  to 

In  (Re  Xi+   J   ImAj.)   <  0 
or  In  (  IXJ   ej^    )   <  0 

or  In  |  Xj     +  J  (  0   ±  2n7T)<0 

which  is  satisfied  if  and  only  if 

IAJ  <  i 

Immediately,  a  direct  analogy  of  the  unit  circle  in  the  X    -  plane  can  be 
drawn  to  the  unit  circle  of  the  z-plane.   This  analogy  can  also  be  demon- 
strated by  considering  the  z- trans  format ion  of  equation  (1-6)  as: 

zX(z)  -  |T  X(z)  +  zX(0) 
(zl  -  $  )X(z)  =  zX(0) 
or,  X(z)  =  z(zl  -  J  )_1  X(0) 

Now,  the  inverse  of  the  matrix  (zl  -J)  is  the  transpose  of  the  cof actor 
matrix,  divided  by  the  determinant  value  of  (zl  -  IE  );  or,  if  the  cofactor 
transpose  matrix  is  denoted  as  [C] 

r*T   xn"1        [C] 
(ZI  "  $  )        ~   |zl  -  £| 


62 


Then 

Since  z  is  iust  a  scalar  variable  over  a  complex  region,  the  determinant 
|  zl   jj  J   is  nothing  more  than  the  characteristic  equation  of  the 
matrix,  whose  characteristic  values  or  eigenvalues  are  the  closed  loop 
root  locations  in  the  z-plane,  or   ^\  -plane.   Similiarly,  in  equation 
(1-1),  the  eigenvalues  of  the  F  matrix  are  the  open  loop  root  locations 
in  the  z-plane  (  \   -plane).   The  vector-matrix  formulation  is  ideally 
suited  for  digital  computation  of  the  state  variable  time  response  in 
discrete  form.   Even  more  important,  it  is  also  suited  for  determination 
of  the  system  dynamic  stability  by  digitally  locating  the  eigenvalues  of 
the  J  matrix  as  a  function  of  any  variable  in  the  2E"  matrix. 
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APPENDIX    II 
FORTRAN    PROGRAM    FOR    MINIMUM    SETTLING    TIME    PARAMETER    PLANE 


•       PROGRAM    DESIGN 
C  PROGRAMMER    L.D.NACF 

C  TH'IS    PROGRAM    COMPUTES    AND    GRAPHS    THE    VALUES    OF    THE 

C  PARAMETERS    KA1    AND    KA2    REQUIRED    FOR    MINIMUM    PROTOTYPE 

C  RESPONSE    OF    ANY    SECOND    ORDER    SAMPLED-DATA    FEEDBACK 

C  CONTROL    SYSTEM    AS    A    FUNCTION    OF    THE    SAMPLING     INTERVAL    T. 

C  INPUT    DATA     IS 

C  CARD    1       VALUES    OF    P,     THE    NEGATIVE    OF    THE    PLANT    S-PLANE 

C  POLE    LOCATION    IN    DECREASING    MAGNITUDE    ORDER»     FORMAT    2F10.5 

C  CARD    2       COEFFICIENTS    OF    THE    S-PLANE    NUMERATOR    OF    PLANT 

C  IN    ASCENDING    ORDER*    FORMAT    3F10.5 

C  CARD    3       T    MIN,    T    MAX*    DELTA    T,     IN    FORMAT    3F10.5.    T    MIN    IS 

C  THE    MINIMUM    VALUE    OF    T    SCANNED.     T    MAX    IS    THE    MAXIMUM 

C  VALUE*    AND    DELTA    T    IS    THE    INCREMENTAL    VALUE    OF    T    USED    IN 

C  CALCULATIONS,    MAXIMUM    NUMBER    OF    POINTS    COMPUTED    IS    500. 

C  T    MIN    MUST    BE    GREATER    THAN    ZERO. 

C  CARD    4       FIRST    LINE    OF    THE    GRAPH    TITLE*    COLUMNS    1-48. 

C  CARD    5       SECOND    LINE    OF    THE    GRAPH    TITLE*    COLUMNS    1-48 

C  CARD    6  EXSCALE»YSCALE    FOR    GRAPH    OUTPUT    IN    FORMAT    2F10.5 

DIMENSION    P(2) »B(3) *ANSK  500) *ANS2(500) »TT(500) 

DIMENSION    ITITLEU2)  .     ^ 

READ    100.P 

100  FORMAT     (2F10.5) 

75  FORMAT     (/////) 

76  FORMAT     (//) 
PRINT    75 
PRINT    101 

101  FORMAT     (47HVALUES    OF    P,    NEGATIVE    OF    S-PLANE    POLE    LOCATIONS) 
PRINT    76 

PRINT    100, (P(I ) ,1=1,2) 
READ    102*8 

102  FORMAT     (3F10.£)  ..) 
PRINT    75 

PRINT     103 

103  FORMAT     (47HPLANT    NUMERATOR    COEFFICIENTS     IN    ASCENDING    ORDER) 
PRINT    76 

PRINT    102*(B( I) ♦ 1=1,3) 
READ    102*    TMIN*TMAX,DELTAT 
PRINT    104 

104  FORMAT     (/////* 3X » 5HT    MIN,6X,5HT    MAX »6X , 7HDELTA    T,//) 
PRINT    102,    TMIN,TMAX,DELTAT 

200  FORMAT     (6A8) 

READ    200,     ( ITITLE( I ) »I=1.6) 
READ    200, ( ITITLE( I )♦ 1=7,12)        v 
PRINT    201 

201  FORMAT     (/////, 12HGRAPH    TITLES,//) 
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PRINT    200,  (  I  TITLE!  I )  »  I  =  l»6) 
PRINT    200, (  ITITLE(  I  >, 1=7,12) 
RFAD    202,    EXSCALE,YSCA|_F 

202  FORMAT  (2F10.5) 
PRINT'203 

203  FORMAT  (/////♦ 10X » 7HX-SCALE * 10X , 7HY-SCALE ♦// ) 
PRINT  2C4,EXSCALE»YSCALE 

204  FORMAT  ( 8X , F10 . 5 , 9X , Fl 0. 5 ) 
LAST  =  0 

P1=P(1) 

P2=P (2) 

BO=B(I) 

B1=B(2) 

B2=6(3) 

K  =  0  .  • 

T=TMIN  .  ; 

105  K=K+1 

IF  (P1+P2)     118,112,118 
118  U=1./EXPF(P1*T) 
V=1./EXPF(P2*T) 
/  PHI l=P2*U/( P2-P1 )+Pl*V/ (P1-P2 )  ' 

PH12=U/(P2-P1 )+V/( P1-P2 ) 
PH2'1  =  (U*V-1.  )*(  (P1*P2)  /(P1+P2  )  ) 
PH22=P1*U/(P1-P2)+P2*V/(P2-P1) 
N  =  0 
IF  (PI)    107, 106*107 

106  N=N+1 

107  IF       (P2)  109,108,109 

108  N  =  N+1 

109  IF     (N-1J  110,111,112 

110  Dl=60/(Pl*P2)+( (B2*P1*P1-B1*P1+B0) / ( P1*P1-P1*P2 ) ) *U 
l+( (6  2*P2*P2-B1*P2+B0)/ (P2*P2-P1*P2 ) ) *V 

D2=( (B1-P1-P2)*(-P1 )+B0-Pl*P2 )*U/(P2-P1 ) 
l+( (Sl-Pl-P2)#(-P2)+80-Pl*P2)*V/(Pl-P2) 
GO    TO    113 

111  Dl=( (B1+B0*T)*P1-B0) /(P1*P1 )+( (B2*P1*P1-B1*P1+30)*U)/(P1*P1) 
D2=B0/Pl-( (B2*P1-B1)*P1+B0)*U/P1 

GO    TO    113 

112  Dl=32+Bl*T+B0*T*T/2# 

D2=B1+80*T  i 

PH11=1. 

PH12=T 

PH21=0. 

PH22=1. 

113  DEN0M=D1*(PH21*D1-PH11*D2)-D2*(PH12*D2-PH22*D1) 
DUMA  1= ( PHI 1  +  PH22 ) * ( PH2 1*D1-PH1 1*D2 ) 

1+D2* (PH11*PH22-PH12*PH21 ) 
DUMA2  =  -D1*(PH11*PH22'-PH12*PH21)'-(PH11  +  PH22)*(PH12*D2-PH22*D1  ) 
ANSI (K)=DUMA1/DEN0M 
ANS2(K)=DUMA2/DENOM 
TT(K)=T 
IF     (TMAX-T)  115,115,114 

114  T=T+DELTAT 
GO    TO    105 
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v        •'  '      ■  !  ,  .vi'K-i^lS-r"   IV  WAV/WOT  Vi-'il'      , 


115 
116 

117 


PR  I 
FOR 
PR  I 
FOR 
CAL 

1  0  ,  0 
CAL 

1 .0  ,'0 
END 
END 


NT  116 

MAT  (/////,10X»1HT*19X»3HKA1,17X»3HKA2»//  ) 

NT  117»   (TT( J) »ANS1 (J) »ANS2( J) »J=1»K) 

MAT(5X,F10.5»10X,F10.5»10X,F10.5) 

L  DRAW  U»TTtANSl»0t0t4HKAl  » I T I TLE » EXSCALE ♦ Y5CALE » 

♦0,0,8,8,1*LAST) 

L  DRAW  (<»TT»ANS2.0tO»4HKA2 

,0,0»8,8,1»LA5T) 


»ITITLE»EXSCALEtYSCALE» 
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